{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.8 Nonlinear minimization problems\n",
    "We consider problems of the form\n",
    "\n",
    "$$\\text{find } u \\in V \\text{ s.t. } E(u) \\leq E(v) \\quad \\forall~  v \\in V.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from netgen import gui\n",
    "from ngsolve import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Scalar minimization problems\n",
    "As a first example we take $V = H^1_0$ and \n",
    "\n",
    "$$\n",
    "E(u) = \\int_{\\Omega} \\vert \\nabla u \\vert^2 + u^4 - fu ~dx.\n",
    "$$\n",
    "\n",
    "The minimization is equivalent to solving the nonlinear PDE:\n",
    "\n",
    "$$\n",
    " - \\Delta u + 4 u^3 = f \\text{ in } \\Omega\n",
    "$$\n",
    "\n",
    "We solve the PDE with a Newton iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from netgen.geom2d import unit_square\n",
    "\n",
    "mesh = Mesh (unit_square.GenerateMesh(maxh=0.2))\n",
    "V = H1(mesh, order=4, dirichlet=[1,2,3,4])\n",
    "u = V.TrialFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To solve the problem we use the `Variation` integrator. Based on the symbolic description of the energy functional, it is able to \n",
    "\n",
    "* evaluate the energy functional (`Energy`)\n",
    "\n",
    "$$E(u) \\qquad (E:V \\to \\mathbb{R})$$\n",
    "\n",
    "* compute the Gateau derivative for a given $u$ (`Apply`):\n",
    "\n",
    "$$A(u)(v) = E'(u)(v) \\qquad (A(u): V \\to \\mathbb{R})$$\n",
    "\n",
    "* compute the second derivative (`AssembleLinearization`)\n",
    "\n",
    "$$(\\delta A)(w)(u,v) \\qquad (\\delta A(w): V\\times V \\to \\mathbb{R})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm (V, symmetric=True)\n",
    "a += Variation ( (grad(u)*grad(u) + u**4-u) * dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Equivalent to:\n",
    "```\n",
    "a += (2 * grad(u) * grad(v) + 4*u*u*u*v - 1 * v)*dx\n",
    "``` \n",
    "(which has the same form as the problems in [the nonlinear example](../unit-3.7-nonlinear/nonlinear.ipynb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We recall the Newton iteration (cf. [unit-3.7](../unit-3.7-nonlinear/nonlinear.ipynb) ) we make the loop:\n",
    "\n",
    "- Given an initial guess $u^0$\n",
    "- loop over $i=0,..$ until convergence:\n",
    "  - Compute linearization: $A u^i + \\delta A(u^i) \\Delta u^{i} = 0 $:\n",
    "    - $f^i = A u^i$ \n",
    "    - $B^i = \\delta A(u^i)$ \n",
    "    - Solve $B^i \\Delta u^i = -f^i$\n",
    "  - Update $u^{i+1} = u^i + \\Delta u^{i}$\n",
    "  - Evaluate stopping criteria\n",
    "- Evaluate $E(u^{i+1})$\n",
    "\n",
    "As a stopping criteria we take $\\langle A u^i,\\Delta u^i \\rangle = \\langle A u^i, A u^i \\rangle_{(B^i)^{-1}}< \\varepsilon$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def SolveNonlinearMinProblem(a,gfu,tol=1e-13,maxits=25):\n",
    "    res = gfu.vec.CreateVector()\n",
    "    du  = gfu.vec.CreateVector()\n",
    "\n",
    "    for it in range(maxits):\n",
    "        print (\"Newton iteration {:3}\".format(it),end=\"\")\n",
    "        print (\"energy = {:16}\".format(a.Energy(gfu.vec)),end=\"\")\n",
    "    \n",
    "        #solve linearized problem:\n",
    "        a.Apply (gfu.vec, res)\n",
    "        a.AssembleLinearization (gfu.vec)\n",
    "        inv = a.mat.Inverse(V.FreeDofs())\n",
    "        du.data = inv * res\n",
    "    \n",
    "        #update iteration\n",
    "        gfu.vec.data -= du\n",
    "\n",
    "        #stopping criteria\n",
    "        stopcritval = sqrt(abs(InnerProduct(du,res)))\n",
    "        print (\"<A u\",it,\", A u\",it,\">_{-1}^0.5 = \", stopcritval)\n",
    "        if stopcritval < tol:\n",
    "            break\n",
    "        Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction (V)\n",
    "gfu.vec[:] = 0\n",
    "Draw(gfu,mesh,\"u\")\n",
    "\n",
    "SolveNonlinearMinProblem(a,gfu)\n",
    "\n",
    "print (\"energy = \", a.Energy(gfu.vec))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again, a Newton for minimization is shipped with NGSolve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve.solvers import *\n",
    "gfu.vec[:] = 0\n",
    "NewtonMinimization(a,gfu)\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Nonlinear elasticity\n",
    "\n",
    "We consider a beam which is fixed on one side and is subject to gravity only. We assume a Neo-Hookean hyperelastic material. The model is a nonlinear minimization problem with \n",
    "\n",
    "$$\n",
    "  E(v) := \\int_{\\Omega} \\frac{\\mu}{2} ( \\operatorname{tr}(F^T F-I)+\\frac{2 \\mu}{\\lambda} \\operatorname{det}(F^T F)^{\\frac{\\lambda}{2\\mu}} - 1) - \\gamma ~ (f,v) ~~ dx\n",
    "$$\n",
    "\n",
    "where $\\mu$ and $\\lambda$ are the Lamé parameters and $F = I + D v$ where $v: \\Omega \\to \\mathbb{R}^2$ is the sought for displacement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import netgen.geom2d as geom2d\n",
    "from ngsolve import *\n",
    "\n",
    "geo = geom2d.SplineGeometry()\n",
    "pnums = [ geo.AddPoint (x,y,maxh=0.01) for x,y in [(0,0), (1,0), (1,0.1), (0,0.1)] ]\n",
    "for p1,p2,bc in [(0,1,\"bot\"), (1,2,\"right\"), (2,3,\"top\"), (3,0,\"left\")]:\n",
    "     geo.Append([\"line\", pnums[p1], pnums[p2]], bc=bc)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# E module and poisson number:\n",
    "E, nu = 210, 0.2\n",
    "# Lamé constants:\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "V = H1(mesh, order=2, dirichlet=\"left\", dim=mesh.dim)\n",
    "u  = V.TrialFunction()\n",
    "\n",
    "#gravity:\n",
    "force = CoefficientFunction( (0,-1) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def Pow(a, b):\n",
    "    return exp (log(a)*b)\n",
    "    \n",
    "def NeoHook (C):\n",
    "    return 0.5 * mu * (Trace(C-I) + 2*mu/lam * Pow(Det(C), -lam/2/mu) - 1)\n",
    "\n",
    "I = Id(mesh.dim)\n",
    "F = I + Grad(u)\n",
    "C = F.trans * F.trans\n",
    "\n",
    "factor = Parameter(1.0)\n",
    "\n",
    "a = BilinearForm(V, symmetric=True)\n",
    "a += Variation(  NeoHook (C).Compile() * dx \n",
    "                -factor * (InnerProduct(force,u) ).Compile() * dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We want to solve the minimization problem for $\\gamma = 5$. Due to the high nonlinearity in the problem, the Newton iteration will not convergence with any initial guess. We approach the case $\\gamma = 5$ by solving problems with $\\gamma = i/10$ for $i=1,..,50$ and taking the solution of the previous problem as an initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(V)\n",
    "gfu.vec[:] = 0\n",
    "\n",
    "Draw (gfu, mesh, \"u\")\n",
    "SetVisualization (deformation=True)\n",
    "\n",
    "res = gfu.vec.CreateVector()\n",
    "du = gfu.vec.CreateVector()\n",
    "\n",
    "for loadstep in range(50):\n",
    "    print (\"loadstep\", loadstep)\n",
    "    factor.Set ((loadstep+1)/10)\n",
    "    SolveNonlinearMinProblem(a,gfu)\n",
    "    Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Allen-Cahn equation\n",
    "\n",
    "The Allen-Cahn equations describe the process of phase separation and is the ($L^2$) gradient-flow equation to the energy\n",
    "$$\n",
    "  E(v) = \\int_{\\Omega} \\varepsilon \\vert \\nabla v \\vert^2~+~v^2(1-v^2) ~ dx\n",
    "$$\n",
    "\n",
    "i.e. the solution to the Allen-Cahn equation solves\n",
    "\n",
    "$$\n",
    "\\partial_t u = \\frac{\\delta E}{\\delta u}\n",
    "$$\n",
    "\n",
    "The quantity $u$ is an indicator for a phase where $-1$ refers to one phase and $1$ to another phase. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The equation has two driving forces:\n",
    "\n",
    "- $u$ is pulled into one of the two minima ($-1$ and $1$) of the nonlinear term $u^2(1-u^2)$ (separation of the phases)\n",
    "- the diffusion term scaled with $\\varepsilon$ enforces a smooth transition between the two phases. $\\varepsilon$ determines the size of the transition layer\n",
    "\n",
    "We use the \"SymbolicEnergy\" feature to formulate the energy minimization problem and combine it with an implicit Euler discretization:\n",
    "\n",
    "$$\n",
    " M u^{n+1} - M u^n = \\Delta t \\underbrace{\\frac{\\delta E}{\\delta u}}_{=:A(u)} (u^{n+1})\n",
    "$$\n",
    "\n",
    "which we can interpreted as a nonlinear minimization problem again with the energy\n",
    "\n",
    "$$\n",
    "  E^{IE}(v) = \\int_{\\Omega} \\frac{\\varepsilon}{2} \\vert \\nabla v \\vert^2~+~v^2(1-v^2) + \\frac{1}{2\\Delta t} \\vert v - u^n \\vert^2 ~ dx\n",
    "$$\n",
    "\n",
    "To solve the nonlinear equation at every time step we again rely on Newton's method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "\n",
    "from netgen.geom2d import *\n",
    "\n",
    "periodic = SplineGeometry()\n",
    "pnts = [ (0,0), (1,0), (1,1), (0,1) ]\n",
    "pnums = [periodic.AppendPoint(*p) for p in pnts]\n",
    "lright = periodic.Append ( [\"line\", pnums[0], pnums[1]],bc=\"periodic\")\n",
    "btop = periodic.Append ( [\"line\", pnums[1], pnums[2]], bc=\"periodic\")\n",
    "periodic.Append ( [\"line\", pnums[3], pnums[2]], leftdomain=0, rightdomain=1, copy=lright, bc=\"periodic\")\n",
    "periodic.Append ( [\"line\", pnums[0], pnums[3]], leftdomain=0, rightdomain=1, copy=btop, bc=\"periodic\")\n",
    "\n",
    "mesh = Mesh (periodic.GenerateMesh(maxh=0.2))\n",
    "V = Periodic(H1(mesh, order=4, dirichlet=[]))\n",
    "u = V.TrialFunction()\n",
    "\n",
    "eps = 4e-3\n",
    "dt = 1e-1\n",
    "gfu = GridFunction(V)\n",
    "gfuold = GridFunction(V)\n",
    "a = BilinearForm (V, symmetric=False)\n",
    "a += Variation( (eps/2*grad(u)*grad(u) + ((1-u*u)*(1-u*u)) \n",
    "                     + 0.5/dt*(u-gfuold)*(u-gfuold)) * dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "gfu = GridFunction(V)\n",
    "gfu.Set(sin(2*pi*x))\n",
    "#gfu.Set(sin(1e8*x)) #<- essentially a random function\n",
    "Draw(gfu,mesh,\"u\")\n",
    "SetVisualization (deformation=True)\n",
    "t = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "for timestep in range(50):\n",
    "    gfuold.vec.data = gfu.vec\n",
    "    SolveNonlinearMinProblem(a,gfu)\n",
    "    Redraw() \n",
    "    t += dt\n",
    "    print(\"t = \", t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Minimal energy extension (postscript in [unit-2.1.3](../unit-2.1.3-bddc/bddc.ipynb) )\n",
    "\n",
    "$$\n",
    "  u \\in V^{ho,disc}, \\quad u^{lo,cont} \\in V^{lo,cont}, \\quad \\lambda \\in V^{lo,disc},\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from netgen.geom2d import unit_square\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
    "fes_ho = Discontinuous(H1(mesh, order=10))\n",
    "fes_lo = H1(mesh, order=1, dirichlet=\".*\")\n",
    "fes_lam = Discontinuous(H1(mesh, order=1))\n",
    "fes = FESpace([fes_ho, fes_lo, fes_lam])\n",
    "uho, ulo, lam = fes.TrialFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "\\int_{\\Omega} \\frac12 \\Vert \\nabla u \\Vert^2  - u + \\sum_T \\sum_{V \\in V(T)} ((u-u^{lo})\\cdot \\lambda)|_{V} \\longrightarrow \\operatorname{min}!\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes)\n",
    "a += Variation(0.5 * grad(uho)*grad(uho)*dx \n",
    "               - 1*uho*dx \n",
    "               + (uho-ulo)*lam*dx(element_vb=BBND))\n",
    "gfu = GridFunction(fes)\n",
    "solvers.Newton(a=a, u=gfu)\n",
    "Draw(gfu.components[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The minimization problem is solved by the solution of the PDE:\n",
    "$$\n",
    "\\int_{\\Omega} \\nabla u \\cdot \\nabla v = \\int_{\\Omega} 1 \\cdot v \\quad \\forall ~ v \\in V^{ho,disc}\n",
    "$$\n",
    "under the constraint\n",
    "$$\n",
    "  u(v) = u^{lo}(v) \\quad \\text{ for all vertices } v \\in V(T) \\text{ for all } T.\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
